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Abstract 

The process of dimerization, in which two monomers bind to each other and form a dimer, is 
common in nature. This process can be modeled using rate equations, from which the average copy 
numbers of the reacting monomers and of the product dimers can then be obtained. However, the 
rate equations apply only when these copy numbers are large. In the limit of small copy numbers 
the system becomes dominated by fluctuations, which are not accounted for by the rate equations. 
In this limit one must use stochastic methods such as direct integration of the master equation 
or Monte Carlo simulations. These methods are computationally intensive and rarely succumb 
to analytical solutions. Here we use the recently introduced moment equations which provide a 
highly simplified stochastic treatment of the dimerization process. Using this approach, we obtain 
an analytical solution for the copy numbers and reaction rates both under steady state conditions 
and in the time-dependent case. We analyze three different dimerization processes: dimerization 
without dissociation, dimerization with dissociation and hetero-dimer formation. To validate the 
results we compare them with the results obtained from the master equation in the stochastic limit 
and with those obtained from the rate equations in the deterministic limit. Potential applications 
of the results in different physical contexts are discussed. 

PACS numbers: 02.50.Fz,02.50.Ey,05.10.-a,82.20.-w 
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I. INTRODUCTION 



Dimerization is a common process in physical, chemical and biological systems. In this 
process, two identical units (monomers) bind to each other and form a dimer (A + A — > A 2 ). 
This is a special case of a more general reaction process (hetero-dimerization) of the form 
A + B — > AB. Dimerization may appear either as an isolated process or incorporated 
in a more complex reaction network. The modeling of dimerization systems is commonly 
done using rate equations, which incorporate the mean-field approximation. These equations 
describe the time evolution of the concentrations of the monomers and the dimers. Assuming 
that the system is spatially homogeneous, these concentrations can be expressed either in 
terms of the copy numbers per unit volume or in terms of the total copy number of each 
molecular species in the system. The rate equations are reliable when the copy numbers of 
the reacting monomers in the system are sufficiently large for the mean-field approximation 
to apply. However, when the copy numbers of the reacting monomers are low, the system 
becomes highly fluctuative, and the rate equations are no longer suitable. Therefore the 
analysis of dimerization processes under conditions of low copy numbers requires the use 
of stochastic methods IlL M. These methods include the direct integration of the master 
e quatl0 n B a, a„ d Monte Ca* (MC, — M fl fl fl g Q. The .aster e quatl0 n 
consists of coupled differential equations for the probabilities of all the possible microscopic 
states of the system. These equations are typically solved by numerical integration. However, 
in some simple cases, the steady state probabilities can be obtained by analytical methods 



10|. The difficulty with the master equation is that it consists of a large number of 



equations, particularly if the dimerization process is a part of a larger network. This severely 
limits its usability for the analysis of complex reaction networks 
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12|. Monte Carlo 



simulations provide a stochastic implementation of the master equation, following the actual 
temporal evolution of a single instance of the system. The mean copy numbers of the 
reactants and the reaction rates are obtained by averaging over an ensemble of such instances. 
In both methods, there is no closed form expression for the time dependence of the copy 
numbers and reaction rates. 

Recently, a new method for the stochastic modeling of reaction networks was developed, 



which is based on moment equations [13| 
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151 ] . The moment equations are much more 



efficient than the master equation. They consist of only one equation for each reactive 
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species, one equation for each reaction rate, and in certain cases one equation for each 
product species. Thus, the number of moment equations that describe a given chemical 
network is comparable to the number of rate equations, which consist of one equation for 
each species. Moreover, unlike the rate equations, the moment equations are linear equations. 
In some cases, this feature enables to obtain an analytical solution for the time dependent 
concentrations. 

In this paper, we apply the moment equations to the analysis of dimerization systems with 
fluctuations. These equations are accurate even in the limit of low copy numbers, where 
fluctuations are large and the rate equations fail. We show how to obtain an analytical 
solution for the time dependent concentrations of the reactant and product species as well 
as for the reaction rate. We identify and characterize the different dynamical regimes of 
the system as a function of the parameters. We examine the validity of our solution by 
comparison to the results obtained from the master equation. The analysis is performed for 
three variants of the system: dimerization, dimerization with dissociation and hetero-dimer 
formation. 

The paper is organized as follows. In Sec. [IT] we analyze the dimerization process using the 
moment equations and provide an analytical solution for the time-dependent concentrations. 
In Sec. II III we extend the analysis to the case in which dimers may dissociate. In Sec. IIVI 
we generalize the analysis to the formation of hetero-dimers. The results are summarized 
and potential applications are discussed in Sec. [V] 

II. DIMERIZATION SYSTEMS 

Consider a system of molecules, denoted by A, which diffuse and react on a surface or 
in a liquid solution. Molecules are produced or added to the system at a rate g (s _1 ), and 
degrade at a rate d\ (s" 1 ). When two molecules encounter each other they may bind and 
form the dimer D. The rate constant for the dimerization process is denoted by a (s" 1 ). 
For simplicity, we assume that the product molecule, D, is non-reactive and undergoes 
degradation at a rate d 2 (s" 1 ). The chemical processes in this system can be described by 
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(1) 



A. Rate Equations 

The dimerization system described above is characterized by the average copy number of 
the monomers, (Na), and by the average copy number of the dimers, (N^). Denoting the 
average dimerization rate by (R) (s _1 ), the rate equations for this system take the form 



<^A = g - M N A )-2(R) 

d AM = _ MNd) + {R) . (2) 

These equations include a term for each process which appears in Eq. (CQ). The factor of 2 
in the reaction term of the first equation accounts for the fact that the dimerization process 
removes two A molecules, producing one dimer. The dimerization rate, (R), is proportional 
to the number of pairs of A molecules in the system, (Na)({Na) — l)/2, where the factor 
of 1/2 is absorbed into the rate constant a. As long as the copy number of A molecules is 
large, it can be approximated by 



(R) = a(N A ) 2 . (3) 



Eqs. (121) form a closed set of two non-linear differential equations. Their steady state solution 
is 



where 7 is the reaction strength parameter given by 



ag 



(5) 




Two limits can be identified. In the limit where 7 > 1, the steady state dimerization 
rate satisfies (R) ss ~ g/2, and the steady state dimer population is (Nd) ss — g/2d2- This 
means that almost all the monomers that are generated end up in dimers and the monomer 
degradation process becomes irrelevant. Therefore, this limit is referred to as the reaction- 
dominated regime. The degradation-dominated limit is obtained when 7 <C 1. In this 
limit (Na) ss — g/d\, (R) ss ~ ag 2 jd 2 = jg and (Nd) ss ~ ag 2 / \d\ 2 ^2) , namely most of 
the monomers that are generated undergo degradation and only a small fraction end up in 
dimers. 

The time-dependent solution for the population size of the A molecules can be obtained 
from the first equation in Eqs. ([2]). The result is 



initial conditions. In the reaction-dominated regime, where ag ^> df, the relaxation time 
converges to r ~ 1/ \/8ag. In the degradation-dominated regime, it approaches r ~ l/d\. 

The rate equation analysis is valid as long as the copy numbers of the reactive molecules 
are sufficiently large In the limit in which the copy number (Na) of the monomers is 
reduced to order unity or less, the rate equations [Eqs. (j2j)] become unsuitable. This limit 
can be reached in two situations: when the monomer concentration is very low or when 
the volume of the system is very small. In the limit of low copy number of the monomers, 
the system becomes dominated by fluctuations which are not accounted for by the rate 
equations. A useful characterization of the system is given by the system-size parameter 



which approximates the copy number of the monomers in case that the dimerization is 
suppressed. The parameter Nq provides an upper limit for the monomer population size 
under steady state conditions. It can be used to characterize the dynamical regime of the 
system. In the limit where iVo ^> 1 the copy number of the monomers is typically large and 
the rate equations are reliable. However, when N < 1 the system may become dominated 



(N A ) = (N A r- — (1 + Ce^) , (6) 



where r 



1 / \f d\ + 8ag is the relaxation time and the parameter C is determined by the 




(7) 
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by fluctuations. In this regime, the rate equations fail to account for the population sizes 
and the dimerization rate. In Fig. [I] we present a schematic illustration of the parameter 
space in terms of 7 and N , identifying the four dynamical regimes. 



B. Moment Equations 

To obtain a more complete description of the dimerization process, which takes the fluc- 
tuations into account, we present the master equation approach. The dynamical variables of 
the master equation are the probabilities P(Na, Nd) of having a population of Na monomers 
and Nd dimers in the system. The master equation for the dimerization system takes the 
form 



dP{N £ ND) = 9 [P(Na-1,N d )-P(N a ,N d )} 

+ d 1 [(NA + l)P(NA + l,N D )-N A P(NA,N D )] (8) 

+ d 2 [(N D + 1)P(N A , N D + 1)- N D P(N A , N d )] 

+ a[(N A + 2)(N A + l)P(N A + 2,N D -1)- N A (N A - l)P(N A , N D )]. 

The first term on the right hand side accounts for the addition or formation of A molecules. 
The second and third terms account for the degradation of A and D molecules, respectively. 
The last term describes the reaction process, in which two A molecules are annihilated and 
one D molecule is formed. The dimerization rate is proportional to the number of pairs of 
A molecules in the system, given by Na{Na — l)/2. Therefore, the dimerization rate can be 
expressed in terms of the moments of P(N A , Nd) as 

(R) = a ((N A 2 ) - (N A )) , (9) 
where the moments are defined by 

00 

(N n A N™) = N2N%P(N A ,N D ), (10) 

N A =0 
N D =0 

and n and m are integers. Note that in the stochastic formulation, the expression used 
for the dimerization rate, (R), is different than in the deterministic approach [Eq. (j3J)]. 
The two expressions are equal in case that P(Na) is a Poisson distribution, for which the 
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mean and the variance are equal. The master equation (jSj) for the dimerization system can 
be analytically solved to obtain the steady state probabilities P(Na)- This solution can be 



found at refs. 



lOj . However, an analytical solution for the time dependent case is currently 



not available. For dimerization systems in the degradation-dominated limit, the probability 
distribution P(Na) approaches the Poisson distribution. However as 7 increases, and the 
system enters the reaction-dominated regime, P(Na) becomes different from Poisson. In 
Fig. [2] we present the marginal probability distributions P(Na) (circles) as obtained from 
the master equation for four choices of the parameters, each in one of the four regimes 
shown in Fig. [U In the limit of ^0 < 1 and 7 ^> 1 (a) the system is in the reaction 
dominated regime (quadrant I in Fig. [TJ, and the results obtained from the master equation 
deviate from Poisson (solid lines). Here the parameters are g = 0.5, d\ — 2, a — 200 and 
d.2 = 10 (s -1 ). The distribution shown in (b), where iVo > 1 and 7 ^> 1 (quadrant II in 
Fig. [1]), also deviates from the Poisson distribution. Here the parameters are g = 100, 
d± = 1, a — 1 and 62 = 10 (s _1 ). In the limit of iVo <C 1 and 7 < 1 (c) the system is 
in the degradation-dominated regime (quadrant III in Fig. [T]), and correspondingly P(Na) 
coincides with the Poisson distribution. Here the parameters are g = 0.5, d\ = 5, a = 1 
and d 2 = 10 (s" 1 ). Finally, in (d) iVo ^> 1 and 7 1 (quadrant IV in Fig. [1]), the system 
is dominated by degradation and as before, the master equation results coincide with the 
Poisson distribution. Here the parameters are g — 10, d± — 1, a — 5 x 10 -3 and d 2 = 10 



In addition to the analytical solution mentioned above, the master equation [Eq. (jSJ)] can 
also be integrated numerically using standard steppers such as the Runge Kutta method 
{17I. In numerical simulations, one has to truncate the master equation in order to keep 
the number of equations finite. This is achieved by setting upper cutoffs iV™ ax and Np ax on 
the numbers of A and D molecules, respectively. This truncation is valid if the probability 
for the number of molecules of each type to exceed the cutoff is vanishingly small. 

The population sizes of the A and D molecules and the dimerization rate are expressed in 
terms of the first moments of P(Na, Njj) and one of its second moments, (Na 2 )- Therefore, 
a closed set of equations for the time derivatives of these first and second moments could 
directly provide all the information needed in order to evaluate the population sizes and the 
dimerization rate is|. Such equations can be obtained from the master equation using the 
identity 
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dWNB) = J- N2N-P(N A ,N D ). (11) 



dt 

N A =0 
N D =0 



Inserting the time-derivative P(N A , Np) according to Eq. (jSJ), one obtains the moment 
equations. The equations for the average copy numbers are 



d{NA) = g-d 1 (N A )-2(R) 



dt 
d(N D ) 



-d 2 (N D ) + (R), (12) 



dt 

while the equation for the dimerization rate is 

= {2ag + 4a 2 )(N A ) + (10a - 2d 1 )(R) - 4a 2 (N A 3 ). (13) 

Eqs. ( [121) have the same form as the rate equations ([2]). However the term for the dimer- 
ization rate, (R), as appears in the moment equations is different from the analogous term 
in the rate equations [Eq. 

Eqs. (fT2"l) . together with Eq. (TlBl) . are a set of coupled differential equations, which are 
linear in terms of the moments. Although we have written the equations only for the relevant 
first and second moments, the right hand side of Eq. (TlBl) includes the third moment for 
which we have no equation. In order to close the set of equations one must express this third 
moment in terms of the first and second moments. Different expressions have been proposed. 
For example, in the context of birth-death processes the relation (N A 3 ) = (N A 2 ) (N A ) was 



used 19|. This choice makes the moment equations nonlinear, which might affect their 
stability. Another common choice is to assume that the third central moment is zero (which 
is exact for symmetric distributions) and use this relation to express the third moment in 
terms of the first and second moments [2^]. Here we use a different approach. We set up the 
closure condition by imposing a highly restrictive cutoff on the master equation. The cutoff 
is set at iV™ ax = 2. This is the minimal cutoff that still enables the dimerization process 
to take place. Under this cutoff, one obtains the following relation between the first three 



moments 



13| 



(N a 3 )=3{N a 2 )-2(Na). (14) 



S 



Using this result, one can bring the moment equations [Eqs. (fl2l) - (fT3l) ] into a closed form: 



d -^ = 9 - dl (N A )-2(R) 
d(N D ) 



dt 

d(R) 
dt 



-d 2 (N D ) + (R) 
2ag(N A ) -2{d 1 + a){R). (15) 



Numerical integration of these equations provides all the required moments, from which the 
population sizes and the dimerization rate are obtained. 

1. Steady State Analysis 

The steady-state solution of the moment equations takes the form 



(N A r - 9{a+di] 



2ag + dia + d{ 



(N D )" 



d 2 {2ag + d^a + d\) 

/m- = ^! — (16) 

2ag + d x a + d{ 

In the limit of very small copy numbers the approximation appearing in Eq. (1141) is valid. 
Thus, in this limit the moment equations provide accurate results, both for the population 
sizes (first moments) and for the dimerization rate (involving a second moment). To evaluate 
the validity of the moment equations in the limit of large copy numbers, we compare Eqs. 
( jl~6l) with the solution of the rate equations (jlj), which is valid in this limit. Consider the 
large-system limit, where N ^> 1 [Eq. fifty] . In this limit, the common term in the denom- 
inators in Eqs. (|T6|) approaches di 2 (2'j + 1), where 7 is the reaction strength parameter, 
given by Eq. ([5]). Thus, in the degradation-dominated limit, where 7 1, the steady state 
solution of the moment equations approaches 
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(N D )" = a- S 



d 2 ^i 

(fl)- = X. (17) 
0(1 

Here we use the fact that in order to satisfy both the large-system limit (Nq 3> 1), and the 
degradation-dominated limit (7 < 1), one must also require d\ ^> a. The results appearing 
in Eqs. ( 1171) are consistent with the results of the rate equations in this limit. We conclude 
that the moment equations are also reliable for large populations under the condition that 
the system is in the degradation-dominated regime. To test the results of the moment 
equations for large systems in the reaction-dominated regime we examine the case of 7 ^> 1. 
Here Eqs. (fTBl) are reduced to 



2a 

(fl>- = I . (18) 

In this limit, the monomer copy number {Na) ss , obtained from the moment equations, does 
not match the rate equation result. Nevertheless, the results for the dimer population size, 
(Nd) ss , and for the dimerization rate (R) ss , do converge to the results obtained from the 
rate equations. 

We thus conclude that the accuracy of the moment equations is maintained well beyond 
the small system limit. The equations provide accurate results for the dimer copy number, 
(Nd) ss , and for the dimerization rate, (R) ss , for both small and large systems. As for the 
monomer copy number, the moment equations provide an accurate description in all limits, 
except for the limit where both iVo > 1 and 7 3> 1 (quadrant II in Fig. [1]). In Tabled 
we present a characterization of the different dynamical regimes and the applicability of the 
moment equations for the evaluation of the copy numbers and the dimerization rate in each 
regime. 

In Fig. [3] we present the monomer copy number (Na) ss (circles), the dimer copy number 
(Nd) ss (squares) and the dimerization rate (R) ss (triangles), as obtained from the moment 
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equations, versus the reaction strength parameter, 7. The rate constants are g = 0.01, d\ = 1 
and d,2 = 5 (s -1 ). The reaction rate, a, is varied. These parameters satisfy the small-system 
limit iVo <C 1. The moment equation results are in excellent agreement with those obtained 
from the master equation (solid lines). However, since the populations are small, the results 
of the rate equations show deviations (dashed lines). In Fig. H]we present (Na) ss (circles), 
(Nd) ss (squares) and (R) ss (triangles), as obtained from the moment equations, versus the 
reaction strength parameter, 7. Here the rate constants are g = 10 3 , d\ — 0.1 and d% — 0.1 
(s _1 ). As before, the reaction rate, a, is varied. These parameters satisfy the large-system 
limit No ^> 1, and thus the rate equation results (dashed lines) are accurate. Although the 
populations are large for the entire parameter range displayed, the results of the moment 
equations are in excellent agreement with those obtained from the rate equations. The only 
deviation appears in the results for the monomer population in the limit 7 ^> 1. For the 
parameters used in this simulation it was impractical to simulate the master equation. 

In any chemical reaction it is important to characterize the extent to which fluctuations 
are significant. From the master equation, one can evaluate the fluctuation level in the 
monomer copy number, given by the variance 

a 2 = (N A 2 ) - (N A ) 2 , (19) 

where a is the standard deviation. The problem is that the expression for a includes the 
first moment (Na), which is not always accurately accounted for by the moment equations. 
However, when the copy number is sufficiently large, the rate equations apply, and thus one 
can extract the value of (Na) 2 in this limit from the rate equations. On the other hand, the 
moment equations account correctly for the second moment (Na 2 ) by (Na 2 ) = (R)/a+ (Na)- 
Using this relation, the result at steady state is 



a 2 



g[d{W{d 1 +g)+a{2d 2 1 +d 1 g+2g 2 
(2ag+adi+dlY 



l53r^-^(-l + v / T + ^) Z for iV >l. 

In Fig. [5]we present the coefficient of variation a j (Na) ss (circles), as obtained from Eq. (1201) 
versus the system size parameter N Q . The parameters are di = 1, a = 1, d 2 = 5 (s" 1 ) and g is 
varied. Here (Na) ss was extracted from the moment equations for N Q < 1, and from the rate 
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for iVo < 1 

(20) 



equations for No > 1. In the small-system limit (Nq <C 1), the average fluctuation becomes 
much larger than (Na) ss - The system is thus dominated by fluctuations. As the system 
size increases, a becomes small with respect to (Na) ss , implying that the system enters the 
deterministic regime. In order to validate our results, we compare them with results obtained 
from the master equation (solid line). For iVo < 1 the agreement is perfect, as in this limit 
the moment equations are expected to be accurate. A slight deviation appears for Nq > 1, 
where a is constructed by combining results obtained from the moment equations and from 
the rate equations. In both limits, Eq. (1201) is found to provide a good approximation for 
the fluctuation level of the system. 

2. Time- Dependent Solution 

The time dependent solution for (Na) can be obtained by solving the two coupled equa- 
tions for (Na) and for (R) in Eqs. (fT5l) . The equation for (Njj) receives input from these 
two equations. However, the dimers are the final products of this network and (./Yd) has no 
effect on (Na) and (R). Thus the equations for (Na) and (R) can be decoupled from the 
equation for (No)- One obtains a set of two coupled linear differential equations of the form 



N = MN + b 



(21) 



where N= ({N A ),{R)), b 



(g, 0) and the matrix M is 




(22) 



The two eigenvectors of the matrix M are 



g: 



iven by 





where u 



\J '4a 2 + d\ + Aad\ — lQag. The corresponding eigenvalues are 



1 




(24) 
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Using the matrix Q = (vi,v 2 ), one can write Eq. (12T1) as 



Q^N = Q 1 MQQ 1 N + Q 1 b. (25) 
The result is a set of two un-coupled differential equations of the form 



/'= I n ° !/"+*, 







T-2 



where / = Q 1 N and k = Q x 6. The solution of Eq. (126]) is 



where Ci and C 2 are arbitrary constants. Multiplying Eq. (j27p from the left hand side by 
the matrix Q one obtains the time dependent solution of Eq. (|2~T1) . which is 



2a.(7 + a«i + d{ 

( R ) = o Z^C I j2 + Q 2 ,iCie"* /Tl + Q 2 , 2 C 2 e-^. (28) 
2ag + adi + af 

The first terms on the right hand side of Eqs. (128]) are the steady state solutions (Na) ss and 
(R) ss as they appear in Eqs. ([16]) . The second and third terms represent the time-dependent 
parts of (Na) and These terms exhibit an exponential decay with two characteristic 
relaxation times, T\ and r 2 . Practically, since T\ < r 2 , the effective relaxation time for the 
monomers is Ta = r 2 . 

In the limit of small copy numbers, where JVo<l, Eqs. fT28l) account correctly for the copy 
numbers and for the reaction rates. In this limit (where g <C di), one obtains ta — \jd\. In 
the limit of large copy numbers, where iVo ^> 1, one has to distinguish between degradation- 
dominated and reaction-dominated systems. Consider a degradation-dominated system, 
where iVo ^> 1 and 7 < 1. These two conditions require that d\ ^> a. As before, the 
effective relaxation time is approximated by ta — This result is consistent with the 

results of the rate equations in this regime [Eq. (jBj)]. A peculiar result arises in the reaction- 
dominated regime, in the limit of large populations. In this limit uj ~ -\/4a 2 — lQag. For 
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a < Ag the parameter u becomes a purely imaginary number. This result leads to spurious 
oscillations in the solution presented in Eqs. fl28|) . As shown for the steady state solution 
[Eq. (1161) ]. the moment equations consistently fail in the limit of reaction-dominated systems 
with large copy numbers. To obtain the relaxation time in this limit, one can rely on the 
results obtained from the rate equations, which give ta ~ 1/y/Sag. The relaxation times in 
all the different limits are summarized in Table [Til 

Finally, we refer to the time evolution of the dimer population (Nd)- The equation for 
(Nd) is the second equation in Eqs. (fT5|) . where (R) is to be taken from Eqs. (|28l) . The 
solution of this equation takes the form 

(N D ) = (N D ) SS + C ie - t/T ' + C 2 e~ t/T2 + C 3 e~ t/T3 , (29) 

where (Nd) ss is taken from Eqs. (EE)), C; = Q 2 ,iCi/[d 2 — (l/r^)], T3 = l/d 2 and C3 is an 
arbitrary constant. The effective relaxation time for the copy number of the dimer product 
depends on the value of T3. If T3 < ta, the copy number of the dimers relaxes rapidly. 
Thus, the time required for the dimers to reach steady state is determined by the monomer 
relaxation time, namely td ~ ta- In the opposite case, where T3 > ta, the monomer 
population reaches steady state quickly, and the production rate of the dimers acts in effect 
as a constant generation rate. Correspondingly, the relaxation time for (Nd) in this limit is 
approximated by td ~ l/d 2 (Table HTj) . 

In Fig. [6](a) we present the time evolution of (Na) (circles) (Nd) (squares) and (R) 
(triangles), as obtained from the moment equations. The parameters are g = 2 x 10~ 3 , 
d\ = 0.05, a = 100 and d 2 = 5 (s -1 ). These parameters correspond to the small system 
limit and to the reaction-dominated regime (quadrant I in Fig. [1]). The moment equations 
(symbols) are in perfect agreement with the master equation (solid lines). The rate equations 
(dashed lines) deviate from the stochastic results both in evaluating the steady state values 
of (Na), (Nd) and (R), and in predicting the relaxation times of (Na) and (R). According 
to the rate equations, this relaxation time should be ta — l/y/8ag ~ 0.8 (s), while according 
to the stochastic description ta — l/d± ~ 20 (s). In Fig. EJ^b) we present results for a system 
in the large population limit and in the regime of reaction-dominated kinetics (quadrant II 
in Fig. []]). Here the parameters are g = 10, d\ = 0.5, a = 1 and d 2 = 10 (s -1 ). Under 
these conditions the moment equations fail to produce the correct time transient, and give 
rise to an oscillatory solution (symbols). In this regime the results obtained from the rate 
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equations (dashed lines) are accurate and coincide with the master equation results (solid 
lines). Note that even in this case the moment equations provide the correct values for 
the dimer production rate, (R), and for the dimer population, (iVo), under steady state 
conditions. In Fig. M^c) the parameters are g = 0.01, d\ — 2, a — 10 and di = 10 (s" 1 ). 
These parameters satisfy the small system limit and are in the kinetic regime dominated 
by degradation (quadrant III in Fig. [1]). The results are in perfect agreement with those 
obtained from the master equation (solid lines). However, the rate equations (dashed lines), 
although displaying similar relaxation times, show significant deviations in the steady state 
values of (Nd) and (R). In Fig. E](d) we present results for the case of a large system, where 
the parameters are g = 10, d\ = 0.5, a = 10~ 3 and di = 0.05 (s" 1 ). These parameters 
correspond to a system in the degradation-dominated regime (quadrant IV in Fig. [[]). 
Although in this system the copy numbers are large, the results obtained from the moment 
equations (symbols) are in perfect agreement with those obtained from the master equation 
(solid lines) and from the rate equations (dashed lines). Here the relaxation time for the 
monomer population is ta — 1/cZi and for the dimer population it is td ~ l/d 2 . 



III. DIMERIZATION-DISSOCIATION SYSTEMS 



To generalize the discussion of the previous Section we now consider the case where the 
dimer product D may undergo dissociation into two monomers, at a rate u (s" 1 ). The 
chemical processes in this system are thus 



A^0 



A + A^D 
D^0 



D ^A + A. (30) 



A. Rate Equations 

The rate equations for this reaction take the form 
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= g - d 1 (N A ) - 2{R) + 2u{N D ) 

= _ { d 2 + u)(N D ) + (R), (31) 

where (R) = a(Nj\) 2 . These equations are similar to Eqs. (J2J), except for the u terms 
which account for the dissociation. We define the effective reaction rate constant as a e ff = 
a[d 2 /{u + d 2 )} such that the effective dimerization rate is {R) c s = a e s(N A ) 2 . Under steady 
state conditions, Eqs. fl3Tj) can be written as 

g-d x {N A ) -2(R) eS = 

(R) cS -d 2 (N D ) = 0. (32) 

They take the same form as Eqs. (j2J), for dimerization without dissociation, under steady 
state conditions. The steady state solution for these equations is 
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d\ 



where 



7cfr = —jr (34) 

is the effective reaction strength parameter. In the limit where d 2 3> u, most dimers undergo 
degradation. The dissociation process is suppressed, and the effective reaction rate constant 
is a e fj ~ a, namely the solution approaches that of dimerization without dissociation. In the 
limit where d 2 u, most of the produced dimers end up dissociating into monomers, and 
correspondingly a c s — >• 0. In this limit, the dimerization and dissociation processes reach a 
balance. The effective dimerization rate vanishes and (Na) ss — g/d±. 

B. Moment Equations 

In order to conduct a stochastic analysis we present the master equation for the 
dimerization-dissociation system, which takes the form 



1(3 



dP(N A ,N D ) 



= g[P(N A - 1, N D ) - P(N A , N D )} 
+ d 1 [(N A + l)P(N A + l,N D )-N A P(N j 



dt 



A,N D )} 



+ d 2 [(N D + l)P(N A ,N D + l)-N D P(N A ,N D )} 



+ a[(N A + 2)(N A + 1)P(N A + 2,N D -1)- N A (N A - 1)P(N A , N D )} 



+ u[(N D + 1)P(N A -2,N D + 1)- N D P(N A , N D )}. 



(35) 



This equation resembles Eq. (jSJ), except for the last term which accounts for the dissociation 
process. The master equation can be solved numerically by imposing suitable cutoffs, iV™ ax 
and Np ax . However an analytical solution is currently unavailable. To obtain a much simpler 
stochastic description of this system we refer to the moment equations. Following the same 
steps as in the previous Section, we impose the minimal cutoffs on the master equation, 
that enable all the required processes to take place. More specifically, we choose iV™ ax = 2 
in order to enable the dimerization. We do not limit the copy number of the dimer, iV£>- 
However, we do not allow N A ^ and No ^ simultaneously, because A and D molecules 
do not react with each other. These cutoffs reproduce the closure condition of Eq. (fl4j) . 
They also gives rise to another closure condition, which is needed here, namely (N a Nd) = 0. 
The closed set of moment equations takes the form 



d{N A ) 



g -d 1 (N A )-2(R)+2u(N D ) 



dt 
d{N D ) 

dt 



(u + d 2 ){N D ) + {R) 



d(R) 



2ag{N A ) - 2(di + a)(R) + 2au(N D ). 



(36) 



dt 



The steady state solution of these equations is 



(Na) 



g(a cS + di) 



(Nd) 



(R) 




(37) 



2ga c s + dia c s + d\ 
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Note that this solution resembles the steady state solution shown in Eqs. ( TT6|) . except for the 
replacement of a by a e g. As before, the validity of the moment equations can be characterized 
by the system size parameter, N , and by the effective reaction strength parameter 7 e g. For 
small systems, where ATqCI, the approximation underlying the moment equations is valid, 
and thus the moment equations provide accurate results for {Na), (No) and (R). In the limit 
of large systems, where Nq ^> 1, the validity of the moment equations can be evaluated by 
comparison with the rate equations. Two limits are observed. In the degradation-dominated 
limit, where 7 e ff -C 1, the solution obtained from the moment equations (I37|) converges to 
the solution obtained from the rate equations ( |33|) . The moment equations are thus valid 
in this limit for the monomer copy number, (Na), as well as for the dimer copy number, 
(Njj), and its production rate, (R). However, for large systems in the reaction-dominated 
limit, where 7 e ff ^> 1, the moment equations converge to the rate equations only for (No) 
and (R). In this limit the monomer population size, (Na), is not correctly accounted for by 
the moment equations. In conclusion, the validity of the moment equations is the same as 
in the case of dimerization without dissociation (Tabled]) under the substitution 7 — » 7 e g. 

In Fig. [7] we present (Na) ss (circles), (Nd) ss (squares) and (R) ss (triangles), as obtained 
from the moment equations for the dimerization-dissociation system versus the effective 
reaction strength, 7 c ff. Here the parameters are g = 0.02, d± — 1, a — 2500 and d,2 = 1 
(s _1 ). The variation of 7 e fj along the horizontal axis was achieved by varying the dissociation 
rate constant, u. For these parameters the system is in the small population limit, namely 
iVo <C 1. The moment equation results are found to be in perfect agreement with the results 
obtained from the master equations (solid lines). However, the rate equations (dashed lines) 
show significant deviations for a wide range of parameters. These deviations are largest 
when the dimerization process is dominant (7 e g > 1) as the effects of stochasticity become 
important. In Fig. [8] we present {Na) 8S (circles), (N D ) SS (squares) and (R) ss (triangles), as 
obtained from the moment equations, versus the effective reaction strength, 7^. Here the 
parameters are g = 1000, d\ — 1, a — 1 and di — 1 (s _1 ). The dissociation rate constant, 
u, was varied. For these parameters the system is in the large population limit, namely 
Nq ^> 1. Although the population sizes of the monomer and of the dimer are large, the 
moment equations are in perfect agreement with the rate equations (dashed lines) in the 
limit of 7 c ff <C 1. For 7 c g ^> 1 this agreement is maintained for the dimer population size 
and for its production rate. In this limit the monomer population size is not accounted for 
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by the moment equations. Slight deviations in (Nd) and (R) appear within a narrow range 
around 7 c g ~ 1. In this narrow range the effective reaction strength parameter is far away 
from either of its limiting values. In any case, these deviations are insignificantly small. 

In the case of the dimerization-dissociation process the moment equations f[3"oT) are a 
set of three linear coupled differential equations. As opposed to the case of dimerization 
without dissociation, here the equation for (Nd) does not only receive input from the other 
two equations, but also generates an output into those equations. This does not enable one 
to solve the first and third equations independently to obtain a time dependent solution 
as shown in the previous Section. Here the time dependent solution will include three 
characteristic time scales for the relaxation times of both (Na), (Nd) and (R). To obtain 
these time scales we first write Eqs. f[3"oT) in matrix form as 



N = MN + b, 
where N = ((N A ), (N D ), (R)), b= (0,0,0) and 



(38) 



/ 



M 



—d\ 2u 

-(u + d 2 ) 
2ag 2au 



y 'lag 'lau —2(d\ + a) j 
The time dependent solution of the moment equations is given by 



(39) 



i=i 



(40) 



where i,j = 1,2,3. Here, N ss = ((Na) ss , (Nd) ss , (R) ss ), and the matrix elements Cy are 
determined by the initial conditions of the system. The relaxation times Tj are 



Ta 



where Xj, j = 1,2,3, are the eigenvalues of the matrix M [Eq. fl39l) ]. The time dependent 
solution obtained from the moment equations applies in the limits where iVo <C 1, or in the 
limits where iVo ^> 1 and 7 c g < 1. In the limit where iVo ^> 1 and 7 c fj ^> 1, the time 
dependent solution should be obtained from the rate equations [Eqs. (T3~Tj) ]. 
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IV. HETERO-DIMER PRODUCTION 



Consider the case where the reacting monomers are from two different types of molecules, 
A and B. Each of these molecules is generated at a rate g A {9 b) and degraded at a rate 
d A {dg). The two molecules react to form the dimer D = AB at a rate a (s^ 1 ). The dimer 
product undergoes degradation at a rate do (s -1 )- For simplicity, here we assume that the 
process of the dimer dissociation is suppressed. The chemical processes in this system are 
thus 






9 A 


.4 





9b 


B 


A 


d,A 





B 







A 


+ B 


a 


D 




0. 



D 

(42) 

The average copy numbers of the reactive monomers and of the dimer product are described 
by the following set of rate equations 



^ = 9A -d A {N A )-(R) 
d(N B ) 



dt 
d(N D ) 



9B-d B (N A ) - (R) 

-d D (N D ) + (R), (43) 



dt 

where (R), the dimer production rate, is given by 

(R) = a(N A )(N B ). (44) 

The master equation for this system describes the time evolution of the probabilities 
P(N A , Nb, Nd) for a population N A molecules of type A, N B molecules of type B and Nd 
dimers D in the system. It takes the form 
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dP{NA ^ B,ND) = g A [P(N A -l : N B ,N D )-P(N A: N B} N D )] 

+ g B [P(N A , N B — 1, Nd) — P(N A , N B , N D )} (45) 

+ d A [(N A + l)P{N A + 1, N B , N D ) - N A P(N A , N B , N D )} 

+ d B [(N B + 1)P(N A , N B + 1, N D ) - N B P(N A , N B , N D )} 

+ d D [(N D + 1)P(N A , N B , N D + 1)- N D P(N A , N B , N D )\ 

+ a[(N A + 1)(N B + 1)P(N A + 1,N B + 1, iV D - 1) - iV A iV B P(iV A , N B , N D )] 

In the stochastic description the production rate of the dimer D is proportional to the 
number of pairs of A and B molecules in the system, namely 

(R) — a(N A N B ). (46) 

A more compact stochastic description can be obtained from the moment equations. Here 
one must include equations for the first moments (N A ), (N B ) and (Nd), and for the pro- 
duction rate, which involves the second moment (N A N B ). The results for the first moments 
can be obtained by tracing over the master equation as shown in Sec. [Til However, when 
deriving the equation for (R) one obtains 



d(R) 
dt 



ag B (N A ) + ag A (N B ) - (d A + d B )(R) 
a 2 ((N A 2 N B ) + (N A N B 2 ) - (N A N B )), 



(47) 



which includes third moments for which we have no equations. To obtain the closure con- 
dition we follow the procedure presented in Sec. [TTJ and impose highly restrictive cutoffs on 
the master equation. Here the cutoff's are chosen as iV™ ax = A r |> Lax = Np ax = 1. These are 
the minimal cutoffs that enable the dimerization process to take place. Under these cutoffs, 
the third order moments appearing in Eq. (1471) can be expressed by 15] 



(N A 2 N B ) = (N A N B 2 ) = (N A N B ). 
One then obtains a closed set of moment equations 



(48) 
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d{N A ) 



9A-d A (N A )-(R) 



dt 



d{N B ) 



g B -d B (N B )-{R) 



dt 



d(N D ) 



d D (N D ) + (R) 



dt 



d(R) 



ag B {N A ) + ag A {N B ) - (d A + d B + a)(R). 



(49) 



dt 



As in the case of the homo-molecular dimerization presented above, the validity of the 
moment equations extends well beyond the cutoff restriction. It can be characterized by 
four parameters. The first two are = g A /d A and Nq = g B /d B , which provide the 
upper limits on the monomer population sizes (N A ) SS and (N B ) SS , respectively. The second 
two parameters are the reaction strength parameters, which in the case of hetero-dimer 
production are 7^ = ag A j \d A d B ) and j B = ag B / \d A d B ). In the limit where the populations 
are small, the moment equations provide accurate results for all the moments appearing in 
Eqs. (1491) . When the populations are large, the moment equations provide accurate results 
for the dimerization rate, (R), and for the dimer population (No)- However, if the reaction 
strength parameters are also large, the moment equations will not correctly account for the 
monomer population sizes, (N A ) and (N B ). 

In Fig. Owe present (N A ) SS (circles), (N B ) SS (squares), (N D ) SS (triangles) and (R) ss (x), 
versus the reaction strength parameters 7,4 = 7^, as obtained from the moment equations. 
Here the parameters are g A = 10 -2 , g B = 10~ 2 , d A = 1, d B = 10, dp = 0.2 (s _1 ), and the 
parameter a is varied. These parameters are within the limit of small populations. The 
results are in perfect agreement with those obtained from the master equation (solid lines). 
The rate equations (dashed lines) show strong deviations, which are mainly expressed in the 
reaction-dominated regime. In Fig. [10] we present (N A ) SS (circles), (N B ) SS (squares), (A^£>) ss 
(triangles) and (R) ss (x), versus the reaction strength parameters 7^ = 75, as obtained 
from the moment equations. Here the parameters are g A = 10 3 , g B = 10 3 , d A = 1, d B = 10, 
do = 0.2 (s _1 ), and the parameter a is varied. These parameters are within the limit of 
large populations. Nevertheless the results obtained from the moment equations for (Nd) ss 
and for (R) ss are in good agreement with those obtained from the rate equations (dashed 
lines) in both the reaction-dominated limit and in the degradation-dominated limit. For the 
monomer population sizes, (N A ) SS and (N B ) SS , the moment equations apply only in the limit 
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where 7 a < 1 and 7^ < 1. 



V. SUMMARY AND DISCUSSION 

We have addressed the problem of dimerization reactions under conditions in which fluc- 
tuations are important. We focused on two types of reactions, ho mo- molecular dimerization 
(A + A — > A 2 ) and hetero-dimer production (A + B — > AB). Common approaches for the 
stochastic simulation of such reaction systems include the direct integration of the master 
equation and Monte Carlo simulations. The master equation involves a large number of cou- 
pled equations, for which there is no analytical solution in the time-dependent case. Monte 
Carlo simulations are often computationally intensive and require averaging over large sets 
of data. As a result, the relaxation times and the steady state populations for given values 
of the rate constants can only be obtained by numerical calculations. 

Here we have utilized the recently proposed moment equations method, in order to obtain 
an analytical solution for the relaxation times and for the steady state populations. The 
moment equations provide an accurate description of dimerization processes in the stochastic 
limit, at the cost of no more than three or four coupled linear differential equations. Another 
useful feature of these equations is that in certain cases they also apply in the deterministic 
limit. Using the moment equations we obtained a complete time dependent solution for 
the monomer population (Na), the dimer population (Nd) and the dimerization rate (R), 
in the case of homo-molecular dimerization. Expressions for the relaxation times and the 
steady state populations were found in terms of the rate constants of the different processes. 
In the case of hetero-dimer production the moment equations include four coupled linear 
equations. These equations can be easily solved by direct numerical integration. However, 
a general algebraic expression for this solution is tedious and was not pursued in this paper. 
Stochastic dimerization processes appear in many natural systems. Below we discuss several 
examples. 



One of the most fundamental chemical reactions taking 



is hydrogen recombination, namely H + H — ► H 2 
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ace in the interstellar medium 
This reaction occurs 



23. 
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on the surfaces of microscopic dust grains in interstellar clouds 25|, |26|, |27(] . The resulting 



H 2 molecu 



es participate in further reactions in the gas phase, giving rise to more complex 



molecules 28J. They also play an important role in cooling processes during gravitational 
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collapse and star formation. In recent years there has been much activity in the computa- 
tional modeling of interstellar chemistry. While the gas phase chemistry can be simulated 
by rate equations l30|, the reactions taking place on the dust grain surfaces often require 
stochastic methodsQ, 4,3]. This is because under the extreme interstellar conditions of low 
gas density, the population sizes of the reacting H atoms on the surfaces of these microscopic 
grains are small and highly fluctuative 



331 ] . The processes taking place on the 



grains are the accretion of H atoms onto the surface, the desorption of H atoms from the sur- 
face, and the diffusion of atoms between adsorption sites on the surface. These processes can 
be described by the dimerization system discussed in Sec. [Ill In recent years, experimental 
work was carried out in an effort to obtain the relevant rate constants and for certain grain 



compositions these constants were found 
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38]. The solution of the moment 



equations, as appears in Sec. [TTl provides the production rate of molecular hydrogen on 
interstellar dust grains, in the limit of small grains and low fluxes, where fluctuations are 
important. 

In the biological context, regulation processes in cells can be described by networks of 



interacting genes [3J 



The interactions between genes include transcriptional regulation 



41|. Due to the small size o 



the eel 
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s, some 
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45]. 



processes as well as protein-protein interactions 

of these proteins my appear in low copy numbers, with large fluctuations 
Deterministic methods are thus not suitable for the modeling of these systems. Dimerization 
of proteins is a common process in living cells. In particular, many of the transcriptional 
regulator proteins bind to their specific promoter sites on the DNA in the form of dimers. 
It turns out that such dimerization, taking place before binding to the DNA, provides an 
effective mechanism for the reduction of fluctuations in the monomer copy numbers 46]. 

In a broader perspective, complex reaction networks appear in a variety of physical con- 
texts. The building blocks of these networks are intra-species interactions and inter-species 
interactions. Thus, the analysis presented in this paper of homo-molecular and hetero- 
molecular dimerization processes, lays the foundations for the analysis of more complex 
networks. Complex stochastic networks are difficult to simulate using standard methods, 
because they require exceedingly long simulation times. The moment equations, applied 
here to dimerization systems, provide a highly efficient method for the simulation of com- 
plex chemical networks. 
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TABLE I: The validity (yQ or invalidity (x) of the moment equations for the evaluation of (Na), 
(Nd) and (R) in the different limits of the dimerization system. The results for (Nd) and (R) 
are valid in all limits. The results for (Na) are invalid in the limit of large systems and reaction- 
dominated kinetics. 
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TABLE II: The relaxation times for (Na) and (Nd) in the different limits of the dimerization 
system. For the monomer population the relaxation time is determined by the degradation rate, 
d\ in three of the four limits. In the reaction-dominated, large system limit the relaxation time is 
determined by y/8ag. In case that the degradation rate of the dimer, c^, is large, its relaxation 
follows that of the monomer population. Otherwise, it is determined by cfo. 
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FIG. 1: Dimerization systems can be classified into four different regimes, characterized by the 
system size parameter, No = g/d\, and by the reaction strength parameter 7 = ag/d\. In the 
large system limit, where Nq S> 1 (quadrants II and IV), the monomer copy number is typically 
large, and the fluctuation level of the system is low. In the small system limit, where iVo <C 1 
(quadrants I and III), the monomer copy number is low, and the system becomes dominated by 
fluctuations. The reaction strength parameter characterizes the dominant dynamical process in the 
system. In the limit where 7 > 1 (quadrants I and II) the process of dimerization is dominant, and 
the degradation is suppressed. In the limit where 7 <C 1 (quadrants III and IV) the degradation 
process is dominant, and only a small fraction of the monomers undergo dimerization. 
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FIG. 2: (color online) The steady state probability distribution P(Na), as obtained from the 
master equation (circles) for (a) small degradation-dominated system; (b) small reaction-dominated 
system; (c) large degradation-dominated system and (d) large reaction-dominated system. For the 
degradation-dominated systems [(c) and (d)], P(Na) can be fitted by a Poisson distribution (solid 
lines). For the reaction-dominated systems [(a) and (b)], P(Na) significantly deviates from the 
Poisson distribution. 
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FIG. 3: (color online) The average monomer copy number, (Na) ss (circles), the average dimer 
copy number, (Njj) ss (squares) and the dimerization rate, {R) ss (triangles), versus the reaction 
strength parameter, 7, as obtained from the moment equations under steady state conditions. The 
parameters used here satisfy the small system limit, namely Nq = 10 -2 . The moment equation 
results are in perfect agreement with the results obtained from the master equation (solid lines). 
However, the rate equation results (dashed lines) show significant deviations. This is because in 
the limit of small monomer copy number the system is dominated by fluctuations, which are not 
accounted for by the rate equations. Note that (Na) ss and (-/Vd) ss are dimensionless while (R) ss is 
units of s" 1 . 
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FIG. 4: (color online) The average monomer copy number, (A^4) ss (circles), the average dimer 
copy number, (Nd) ss (squares) and the dimerization rate, (R) ss (triangles) versus the reaction 
strength parameter, 7, as obtained from the moment equations under steady state conditions. 
The parameters used here satisfy the large system limit, namely Nq = 10 4 . In this limit the rate 
equations (dashed lines) are reliable. The moment equation results are in agreement with the 
results obtained from the rate equations for the dimer copy number and its production rate. This 
is in spite of the fact that the copy numbers are large, far beyond the conditions for which the 
moment equations are designed. For the monomer copy number, the moment equations deviate 
in the reaction-dominated limit (7 > 1). Slight deviations in the results for {Nd) ss and (R) ss 
appear around 7 — 1, where the crossover from the degradation-dominated regime to the reaction- 
dominated regime takes place. 
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FIG. 5: (color online) The coefficient of variation of the monomer copy number, ct/{Na) ss , versus 
the system size parameter Nq (circles), as obtained from Eq. (|20p . As the average population size 
increases, a becomes smaller than {Na) ss . For Nq < 1 the results are in perfect agreement with 
those obtained from the master equation (solid line). For iVo > 1, where the expression in Eq. (|20p 
combines results obtained from the moment equations and the rate equations, a slight deviation 
emerges. Nevertheless, Eq. (|2Up is shown to provide a good approximation for a. 
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FIG. 6: (color online) The time dependence of (Na) (circles), (Njj) (squares) and (R) (triangles), 
as obtained from the moment equations. Four different limits are observed. In the limit where 
iVo <C 1 and 7 > 1 (a) the moment equations are in perfect agreement with the master equation. 
In this limit, the rate equations fail to account both for the steady state copy numbers and for 
the relaxation time. According to the rate equations, the relaxation time is r ~ l/y/8ag, while 
according to the moment equations it is r ~ 1/cZi, in agreement with the master equation. In 
the limit where Nq 3> 1 and 7 ^> 1 (b), the rate equations and the master equation are in good 
agreement. In this limit the moment equations fail to produce the correct time dependent solution, 
as they predict an oscillatory convergence towards steady state. Note, however, that even in this 
limit, the moment equations do provide correct results for the steady state values of {Njj) and 
(R). In the limit where iVo <C 1 and 7 <C 1 (c), the moment equations are in perfect agreement 
with the master equation (solid lines). The rate equations (dashed lines), which do not account 
for fluctuations, fail in this limit. The limit where iVo S> 1 and 7 -C 1 is shown in (d). Although 
the copy numbers are large in this limit, the moment equations are still in perfect agreement with 
the master equation. Not surprisingly, the rate equations also apply. 
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FIG. 7: (color online) The average monomer copy number, {Na) ss (circles), the average dimer copy 
number, (iV*£>) ss (squares) and the dimerization rate, {R) ss (triangles), versus the effective reaction 
strength parameter, 7 e ff , as obtained from the moment equations for the dimerization-dissociation 
reaction. The parameters used here represent the small system limit, namely Nq = 0.02. The 
moment equation results are in perfect agreement with those obtained from the master equation 
(solid lines). However the results of the rate equations (dashed lines) show significant deviations. 
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FIG. 8: (color online) The average monomer copy number, {Na) ss (circles), the average dimer copy 
number, {Nd) ss (squares) and the dimerization rate, {R) ss (triangles), versus the effective reaction 
strength parameter, 7 e ff , as obtained from the moment equations for the dimerization-dissociation 
reaction. The parameters used here represent the large system limit, namely Nq = 10 3 . In this 
limit the rate equations (dashed lines) are reliable. The moment equation results are in agreement 
with those obtained from the rate equations for the dimer copy number and its production rate. 
However, for the monomer copy number, the moment equations deviate in the reaction-dominated 
limit (7 c g > 1). Slight deviations in the results for {Nu) ss and (i?) ss appear around 7 e g ~ 1, where 
the crossover from the degradation-dominated regime to the reaction-dominated regime takes place. 
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FIG. 9: (color online) The average monomer copy numbers, {Na) ss (circles) and {Nb) ss (squares), 
the average dimer copy number, {Nd) ss (triangles), and the dimerization rate, {R) ss (x), versus 
the reaction strength parameters, 7^4 = jb, as obtained from the moment equations, for the 
hetero-dimer production system. The parameters used here represent the small system limit. The 
moment equation results are in perfect agreement with those obtained from the master equation 
(solid lines). However the results of the rate equations (dashed lines) show significant deviations. 
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FIG. 10: (color online) The average monomer copy numbers, {Na) ss (circles) and {Nb) ss (squares), 
the average dimer copy number, (Nd) ss (triangles), and the dimerization rate, {R) ss (x), versus 
the reaction strength parameters, -ja = 7_b 5 as obtained from the moment equations for the hetero- 
dimer production system. The parameters used here represent the large system limit. In this 
limit the rate equations (dashed lines) are reliable. The moment equation results are in agreement 
with the results obtained from the rate equations for the dimer population and its production rate. 
However, for the monomer copy numbers, the moment equations deviate in the reaction-dominated 
limit (j A (B) > !)■ 
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